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Abstract. The quantum dissipative dynamics of a tunneling process through double 
t^ ' barrier structures is investigated on the basis of a rigorous treatment. We employ a 

Caldeira-Leggett Hamiltonian with an effective potential calculated self-consistently, 

accounting for the electron distribution. With this Hamiltonian, we use the reduced 
' ^ ■ hierarchy equations of motion in the Wigner space representation to study non- 

Markovian and non-perturbative thermal effects at finite temperature in a rigorous 
CJ ' manner. We study the current-voltage relation of the resonant tunneling diode for 

several widths of the contact region, which consists of doped GaAs. Hysteresis and 

both single and double plateau-like behavior are observed in the negative differential 
K* ' resistance (NDR) region. While all of the current oscillations decay in time in the 

NDR region in the case of a strong system-bath coupling, there exist non-transient 
\/^ ' oscillations in some parts of the plateau in the NDR region in the case of weak coupling. 

VO ■ We find that the effective potential in the oscillating case possesses a basin-like form on 

^^ I the emitter side (emitter basin) and that the current oscillation results from tunneling 

C3 ■ between the emitter basin and the quantum well in the barriers. We find two distinct 

types of current oscillations, with large and small oscillation amplitudes, respectively. 

The results of eigenstates analysis indicate that the first type is caused by a transition 

between ground tunneling states and adjacent excited states in the emitter basin, while 
j^ , the second type is caused by a transition between intermediate tunneling states and 

higher states. These two types of oscillation also appear differently in the Wigner 

space, with one exhibiting tornado-like motion and the other exhibiting a two piston 

engine-like motion. 
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1. Introduction 

Quantum coherence and its destruction by coupling to a dissipative environment play 
an important role in the transport phenomena of a particle moving in a potential. Well 
known examples include electron transfer in molecular and biological systems [U El [U H] , 
many chemical reactions [SI El [TJ [H], SQUID rings [HI [10], quantum ratchets [HI fT2] . 
nonlinear optical processes [131 [E]) cind tunneling processes in device systems [151 US]- 
Such systems are commonly modeled as one- dimensional or two-dimensional potential 
systems coupled to heat bath degrees of freedom, which drive the systems toward 
the thermal equilibrium state. The heat bath degrees of freedom are then reduced 
using such methods as the projection operator method or the path integral method, for 
example. Many equations of motion have been derived for the purpose of understanding 
the quantum aspects of dissipative dynamics [T71[T8l[T9l[2nl[2Tl[22l[28l[2il[25]. 
Because a complete picture of quantum dissipative dynamics must treat phenomena 
that can only be described in real time, a great deal of effort has been dedicated 
to the problem of numerically integrating these equations of motion in real time 
[SSlEllEllESlEniETlEHlEg]. Although these equations are analogous to the classical 
kinetic equations, which have proved to be useful for classical transport problems, such 
equations cannot be derived in a quantum mechanical framework without significant 
approximations and/or assumptions. For example, the quantum Boltzmann equation is 
based on the assumption that the effects of collisions between electrons can be described 
by the rates determined from Fermi's golden rule, and hence it is regarded as a semi- 
classical equation [T71 [HI [19] . Similarly, the quantum Fokker-Planck equation can be 
derived from the Caldeira-Leggett Hamiltonian under a Markovian approximation, but 
in order for this to be possible, the heat bath must be at a sufficiently high temperature, 
in which case most of the important quantum dynamical effects play a minor role 
[201 |2T] . Thus, treatments of the above kinds cannot be used to construct fully quantum 
mechanical descriptions of broad validity. In this paper, we present a rigorous quantum 
mechanical treatment that solves this problem. This treatment employs the reduced 
hierarchy equations of motion (HEOM), and it can be used in application to systems 
for which fully quantum mechanical description is necessary. In particular, the reduced 
HEOM approach can be used to numerically treat non-Markovian system-bath coupling 
in a non-perturbative manner. Here, we apply this approach to study the dynamics 
of a resonant tunneling diode (RTD) described by the Caldeira-Leggett Hamiltonian 
[30l[3ll[32l[33l[3ll[35l[36l[371[38l[39l[10]. 

The RTD system that we consider is modeled by a double barrier structure with an 
electrostatic potential representing a region consisting of an undoped layer positioned 
between two doped layers. A second potential, arising from the charge distribution, 
is determined self-consistently [151 [IB] • Due to quantum effects, RTD systems exhibit 
novel negative differential resistance (NDR) in the current- volt age (I-V) relation [JT]. 
Moreover, current oscillations [42], plateau-like behavior and hysteresis of the I-V curve 
have been observed in the NDR region ^3]. Goldman et al. attributed these phenomena 
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to intrinsic charge bistability. However, Sollner claimed that they result merely from 
resonance with the external circuit [H]. Resonant tunneling diodes are presently the 
highest-frequency active semiconductor devices in existence [l5lll6lll71llHlll9]. Although 
oscillation frequencies above 1 THz have recently been realized [501 |5T1 [52] , the causes of 
such phenomena have not yet been definitely identified. In these phenomena, dissipation 
and fluctuations arise from interaction with the heat bath, which plays an essential role 
as the origin of resistance. 

Frensley et al. discovered NDR in the I-V curve through a numerical computation 
treating a quantum Liouville equation in the Wigner representation that ignored 
phonon-scattering processes [53l Ell EHl EG] [571 EH]- Kluksdahl et al. incorporated 
dissipative and self-consistent effects, employing the Poisson-Boltzmann equation 
by adopting a relaxation time approximation, and succeeded in modeling the 
experimentally observed hysteresis behavior of the I-V curve [60]. Jensen and Buot 
developed a numerical scheme to treat systems of the same kind and found evidence 
that the current oscillation and plateau-like behavior arise from intrinsic bistability 
[6T| [62] . However, because there exists no well-established methodology that can be 
applied rigorously to this type of model and includes the effect of dissipation, which 
is the origin of Joule heat, its application in describing the dynamics of systems under 
such conditions has been quite limited. 

In a previous study [63] , we demonstrated through consideration of a RTD that with 
the HEOM approach, the effects of non-Markovian thermal fluctuations and dissipation 
at finite temperature on quantum transport can be investigated rigorously. We found 
that while most of the current oscillations decay in time in the NDR region, there exists 
persistent oscillation in a plateau of the NDR region. In this paper, we explore the 
cause of current oscillation by considering several widths of the doped contact regions 
and two different cases for the strength of the system-bath coupling. The value of the 
noise correlation time is also chosen to be twice shorter than the previous study. 

This paper is organized as follows. In Sec. 2, we introduce the reduced hierarchy 
equations of motion applicable to the resonant tunneling problem. We then present the 
computational details for the numerical simulations in Sec. 3. Numerical results for the 
I-V curves and current oscillations are presented in Sees. 4. Section 5 is devoted to 
concluding remarks. 

2. Formulation 

We consider the following Caldeira-Leggett Hamiltonian p], which describes the 
dynamics of an electron subjected to a thermal environment: 

^2 
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Here, m, p, and q are the mass, momentum and position variables of the electron, and 
mj,pj,Xj and Uj are the mass, momentum, position and frequency variables of the jth 
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phonon oscillator mode. In equation ([T]), the electron-phonon interaction is given by 

Hi = -Viq)J2a,x,. (2) 

j 

Here, V{q) is any function of q and the quantities aj are coefficients that depend on the 
nature of the electron-phonon coupling. 

The heat bath can be characterized by the spectral distribution function, defined 
by 

and the inverse temperature, /3 = I/ZcbT, where k-Q is the Boltzman constant. We 
assume the Drude distribution, given by 

where the constant 7 represents the width of the spectral distribution of the collective 
phonon modes and is the reciprocal of the correlation time of the noise induced by 
phonons. The parameter ( is related to the electron-phonon coupling strength. For the 
collective heat bath coordinate X = J2j (^j^j-, the canonical and symmetrized correlation 
functions, respectively defined by ^(t) = /3(X;X(t))B and C{t) = i(X(t)X(0) + 
X{0)X{t))B, where X{t) is the Heisenberg representation of X, and (■ ■ ■)b represents 
the thermal average over the bath modes, are given by [311 [30] 



and 






m = -^e-'*'. (5) 



C(t) = Coe-^l*l+5:c,e-'^'=l*'. (6) 

fc=i 



Here, Uk = 27[k/l3h are the Matsubara frequencies, and we have 



Co 
and 



2 ^ 4/3^7 

/3^7 ti (/3M' - (27rA:) 



(7) 



^''~ 2 (/3%)2-(27rfc)2- ^^ 

The function C{t) is analogous to the classical correlation function of X{t) and 
corresponds to the correlation function of the bath- induced noise, whereas \E'(t) 
corresponds to dissipation. The noise C(t) is related to \E'(t) through the quantum 
version of the fiuctuation-dissipation theorem, C[u] = hu coth{f3hco/2)/2'^[u], which 
insures that the system exists in the thermal equilibrium state for finite temperatures 
|64j . Note that in the high temperature limit, Ph'j <^ 1, the noise correlation function 
reduces to C(t) oc e"'^'*'. This indicates that the heat bath oscillators interact with the 
system in the form of Gaussian-Markovian noise. 
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To derive the equation of motion for the electron, we use the reduced density 
operator of the system by taking the trace over the heat bath degrees of freedom: 

p{t) = TrBPtot(t). (9) 

In the path integral representation, the reduced density matrix elements are written 
p{q, g'; t) = J D[q{r)] J D[q'{T)] J dq, j dq[p{qi,q[)pcs{q^ q\ t; qi,q'i} 

where S^lq; t] is the action for the system's Hamiltonian Ha = p'^/'2m + f/(g; t), 

■1 
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p(qi, q'i) is the initial state of the system at time tj, F[q, q'; t] is the influence functional 
[65] . and pcs{q, q', t] qt, q'i) is the initial correlation function between the system and 
the heat bath [66]. The functional integrals for g(r) and q'{T) are carried out from 
qi'ti) = qi to q{t) = q and from q'(ti) = q^ to q'(t) = q', respectively. In our approach, 
we can set pcs{q, q', i] qi) q'i) = 1 and regard p(gj, gQ as the initial condition without 
loss of generality. This is due to the fact that the effects of the initial correlation can 
be taken into account through implementation of a hierarchy of initial conditions that 
can be evaluated numerically. The influence functional for the inverse temperature /3 is 
given by [651 EE] 

F[q, q'; t] = exp U-^^' £ drV^q, q'; r) 



X 



rdT''-^{r-r')V%q, q'; r') + f dr'C{r - r')V\q, g'; r') 



-dr 

(12) 

where V (g, q'; r) = V (g(r)) - V (g'(r)) and V°{q, q'; r) = V (g(r)) + V (g'(r)). If 
we choose K so as to satisfy uk = 2ttK/{[5K) ^ Wc, where Uc is the characteristic 
frequency of the system, the factor e"*^*'*' in equation ([6]) can be replaced with Dirac's 
delta function, using the approximation i/^e"^''!*! ~ 5{t) (for k > K + 1). Therefore, 
C{t) can be expressed as 

C(t) = coe-^l*l + ^c,e-''^l*l+5(t) f] ^. (13) 

fe=l k=K+l ^k 

The reduced hierarchy equations of motion (HEOM) can be obtained by considering 
the time derivative of the reduced density matrix with the kernel given in equations ([5]) 
and (fT3l) . The HEOM have been used to study chemical reactions [321 |33l [671 EH], 
linear and nonlinear spectroscopy [3l|35l|36ll371[3Hll39lll0l|69l[70l[n], exciton 
transfer [72l [73l EH ES |76l [77], electron transfer [711 [791 [801 [H], quantum dots 
[821 [83] . and quantum information [HH [851 US]. A variety of numerical techniques have 
been developed for the HEOM approach in order to accelerate numerical calculations 
[871 [881 [891 [901 [SH [92]. 
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The HEOM are ideal for studying quantum transport systems, in conjunction with 
the Wigner representation, characterized by the Wigner distribution function. 



Wip,q;t) 



_iPi f r r 



(14) 



because they allow us to treat continuous systems utilizing open boundary conditions 
and periodic boundary conditions [551 ES]. The Wigner distribution function is the 
quantum analogue of the classical distribution function in the phase space [531 Ell ESI 
56l [571 [58l l59] . Its classical limit can be computed readily. This is helpful, because 
knowing the classical limit allows us to identify the purely quantum mechanical effects 
[3211331110]. 

Although we can handle any form of V{q) [351 EEl ETl EH [391 HO], here we 
consider the liner-liner system bath coupling case defined by V{q) = q. In the Wigner 
representation, the equations of motion are expressed in hierarchical form as follows 
MM- 

K 
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(15) 



for non- negative integers n,ji,... ,Jk, where we have chosen K such that Uk 3> oJc- In 
equation (ITSl) . — I/qm is the quantum Liouvillian in the Wigner representation, given by 

1 r'^ dp' 



L,^W{p,q)^-^^W{p,q) 



h J-OD 2TTh 



with 



Uw{p,q]t) =2 J dr sin { — 



U(q+-,t 



U^{p-p',q;t)W{p',q)ilQ) 



(17) 



U(q--,t 



The other operators appearing in equation (TT5J1 are the bath-induced relaxation 
operators, defined as 

*-!;. (18) 



SosC 



V H cot -— 

^ 2 \ 2 dp 



and 









/3^7 



cot 



(^)-S4 



(19) 
(20) 

(21) 



Self- excited current oscillations in a resonant tunneling diode 

{3} 

^^ Dontacts I I spacer [ayera ^B barriers (AIGsAs} | | quantunn well 

:■« doped >« Ljndoped >H doped >i 



^^^^^n+S^^^^ 


■ ■ 





(b) 



0.3 


1 












-0.2 
^ 0.1 

LU 

0.0^ 










] 10 20 




30 


40 



Positbn (nrm) 

Figure 1. (a) The structure of the RTD. The weh consists of undoped GaAs (4.520 
nm), the barriers consist of undoped AlGaAs (2.825 nni), the spacer layers consist 
of undoped GaAs (2.825 nm), and the contact regions consist of doped GaAs (with 
a doping concentration of 2 x 10^® cm~'^). In order to elucidate the dependence of 
the current on the width of the contact regions, we carried out computations for three 
values of this widths: (i) 16.950 nm; (ii) 33.900 nm; (iii) 42.375 nm. This figure depicts 
the case of 16.950 nm. (b) The structure of the conduction band edge. The height of 
the potential barriers is 0.27 eV. 



In the case that the quantity N = n + Y.k=iJk satisfies the relation N ^ 
ujc/in.m{'y,l/f3h), this finite hierarchy can be truncated with neghgible error by the 
terminator EOl 



d 



r{n) 



(n) 



^w^Z...(t) = -iL^ + ^W^!..,M 



dt 



(22) 



.(0) 



Note that only Wq'q{p^ q; t) = W{p, q; t) has physical meaning, and the other elements 
{p,q;t) with (n; ji, . . . , j^) ^ (0;0,...,0) are auxihary operators introduced 



W, 



(n) 



J1,---,Jk 



to avoid the explicit treatment of the inherent memory effects that arise in the time 
evolution of the reduced density matrix. In the classical limit, /i — )■ 0, the Wigner 
distribution function corresponds to the classical distribution function. If the noise 
correlation is very short (7 — ?■ 00) and the temperature is high (i.e., the noise is white), 
the quantum Fokker-Planck equation can be derived in a form similar to that of Kramers 
equation [201 EI]- In the present case, however, we cannot employ the white noise 
approximation, because quantum effects play a dominant role in the low temperature 
regime {Phuc 2> 1). Although the original quantum Fokker-Planck equation, derived 
by Caldeira and Leggett, is valid only for sufficiently high temperatures {Phuc ^ 1) 
[20] . our formulation is valid for arbitrary temperature, and it can be applied without 
approximation. 

In equation (TT71) . U{q;t) is the effective potential for the electron, which can be 
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written [SOI ED E2] 

U{q; t) = Ustaticiq) + U^eifiq; t), (23) 

where Ustatic{,q) and Useif{q',t) are the static and self-consistent parts, respectively. As 
the static potential, we employ the double-barrier structure depicted in figure 1. The 
self-consistent part, Useif{q',t) = —e(j){q;t), is calculated from the electron distribution 
at each time step in the integration of equation ( ITSl) using the Poisson equation 



-^[e0(g;t)] = e[n+(g)-P(g,t)], (24) 

where e is the dielectric constant, n+(g) is the doping density and P{q;t) = 
/f^ (ip/(27r^)Wo . IoIP; '?j ''') is the electron density calculated from the Wigner 
distribution. Coupling the HEOM to the Poisson equation, we obtain a fully self- 
consistent model of quantum electron transport. This allows us to examine charge 
redistribution effects. 

3. Computational details 

The equations of motion given in (TTSll were numerically evaluated using finite mesh 
representations of the Wigner distribution functions. The spatial derivative of the 
kinetic term in the Liouville operator, —{p/m)dW{p, q)/dq, was approximated by using 
a third-order left-handed or right-handed difference scheme, depending on the sign of 
the momentum, given by 

^5^ = g^a^fe.*«)+3;yfe,,,) 



- 6W(j,t,qi-i) + W{pt,qj-2)) (25) 



for pk > 0, and 

dW{pk,qj) 



W {pk, qj+2) + 6iy (pfc, qj+i) 



dq 6Ag 

- m{pk,qj) + 2W{pk,q,^r)) (26) 

for Pk < 0, in order to facilitate implementation of the inflow and outflow boundary 
conditions [531 EU [55]. Note that the first-order difference scheme is not sufficiently 
accurate for the present problem, because this scheme introduces artificial viscosity for 
the wavepacket dynamics, which damps the self-excited current oscillations. 

The infiow boundary conditions were set by stipulating that Wj^ ..^^^{p < 0,q = L) 
and Wji,..j^.iP > 0,g = 0) are given by the equilibrium distribution of a free particle 
calculated from the HEOM with periodic boundary conditions. Due to fiuctuations and 
dissipation, the fiow of a wavepacket reaches a steady state even when there exists a non- 
zero bias voltage. The validity of the boundary conditions was verified by considering 
several system sizes. 
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Other derivatives with respect to p were approximated using the fourth-order 
centered difference scheme given by 

- 8W{pk.i,q,) + W{pk-2,qj)), (27) 



and 



d'W{pk,q,] 



W {pk+2, qj) + lQW{pk+i, qj) - 30W{pk, qj) 



dp'^ 12Ap2 

+ l6W{pk-u qj) - W{pk-2, q,)) . (28) 

The mesh size for the position, Ag, and momentum, Ap, are respectively 0.2825 nm 
and 3.540 nm~^. 

We set the parameters used in the HEOM as 7 = 24.2 THz (7-^ = 4.13 fs), C = 72.5 
GHz (C^^ = 13.8 ps), and T = 300 K in order to create conditions close to those used in 
previous theoretical studies. In Appendix B, we also report the results of calculations 
for the strong coupling case, with C, = 120.8 GHz {(~^ = 8.28 ps), to elucidate the role of 
dissipation. The depth of the hierarchy and the number of Matsubara frequencies were 
chosen as A^ G {2 — 6} and K e {1 — 3}, respectively. The effective mass of the electron 
was assumed to be constant across the device and equal to 0.067mo, where mo is the 
electron mass in vacuum. The dielectric constant in equation (^^ was set as e = 12.85. 

As the static double-barrier potential, which models the hetero-structure of GaAs 
sandwiched between two thin AlGaAs layers, we set the widths of quantum well 
(undoped GaAs), barrier (undoped AlGaAs), and spacer layer (undoped GaAs) to be 
4.520 nm, 2.825 nm, and 2.825 nm, respectively. The height of the potential barriers 
was 0.27 eV. The conduction band edge consists of a single quantum well bounded by 
tunneling barriers (figure 1 (b)). The widths of the contact regions (the yellow parts in 
figure 1(a), where GaAs is doped with a concentration of 2 x 10^^ cm~^) were chosen as 
16.950 nm, 33.900 nm, and 42.375 nm. 

In a previous study [63], we chose a smaller value of 7, 7 = 12.1 THz (7^^ = 8.26 
fs), and fixed the width of the contact region as 42.375 nm. In that case, we found 
hysteresis, double plateau-like behavior, and self-excited current oscillation in the 
negative differential resistance (NDR) region of the current-voltage curve. We found 
that while most of the current oscillations decay in time in the NDR region, there exists 
a non-transient oscillation characterized by a tornado-like rotation in the Wigner space 
in the upper plateau of the NDR region. In this paper, we explore the cause of such 
current oscillations by considering several values of the width of the contact regions in 
cases of both weak and strong coupling, characterized by different values of C- 

4. Results 

We determined the characteristics of the current- voltage (I-V) according to the following 
procedure. First, we integrated equation (fT5|) at zero bias voltage without the self- 
consistent part of the effective potential under the inflow boundary conditions specified 
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above. When we obtained the temporal steady state, the obtained distribution was 
then used as the initial distribution for the self-consistent calculation, and we then 
integrated equation ( TT5|) again, with the effective potential t/(g; t) evaluated iteratively 
using the Poisson equation given in i^^. Under this procedure, when the distributions 
reached the genuine steady state Wj^ j^{p,q;t — )■ cxd), the current was calculated by 
I{t) = J dppWQ^,,,fi{p,q;t)/27!-m'^h and then the genuine state was used as the initial 
distributions for the next bias step. Note that if we were to omit the first step of 
the calculation, which does not contain the self-consistency procedure, there would be 
cases in which the distribution does not reach the steady state, in particular when 
the contact regions are wide. Following the above steps, we increased the bias from 
O.OOOV to 0.500V, and then decreased it to O.OOOV with bias steps of 0.01 V in the 
normal region and 0.002 V in the negative differential resistance (NDR) region. At each 
step, we integrated the equation of motion until the system reached the steady state 
distribution. However, in some cases, in the NDR region, current oscillation arose, and 
a steady state was thus not realized. In such cases, the value of current was evaluated 
as a time average between 30 ps to 40 ps. 

In figure 2, we present (a) current- voltage (I-V) relations and (b) the time evolution 
of the self-excited current oscillation in the weak coupling case {( = 72.5 GHz) for 
three values of the width of the contact regions (the yellow parts in figure 1(a)). In 
these plots, the sizes of the contact regions are (i) 16.950 nm, (ii) 33.900 nm, and (iii) 
42.375 nm, respectively, with a fixed doping concentration of GaAs. These graphs reveal 
NDR behavior, hysteresis, and plateaus in the I-V curve. Moreover, self-excited current 
oscillation appears in some regions of the plateau. While Jensen and Buot observed 
only a single plateau similar to that in figure 2 (i-a) with current oscillation [61], our 
results in figure 2(ii-a) and 2(iii-a) exhibit a double plateau structure. In the case of a 
single plateau structure, the width of the plateau is large for a narrow contact region, 
while in the case of a double plateau structure, the width is small for a wide contact 
region. For the sake of comparison, the I-V curves obtained in the strong coupling case 
{( = 120.8 GHz) are presented in Appendix B. In that case, NDR behavior, hysteresis, 
and a single plateau in the I-V curve are observed, but steady current oscillation does 
not appear, even in the plateau. 

As in the case considered in the previous paper [63], we find that most of the 
current oscillations decay in time in the NDR region, but there also exist non-decaying 
oscillations in some regions of the plateau, as seen in figures 2 (i-b) - (iii-b). The Fourier 
components of each persistent oscillation are plotted in figures 3 (i-a) - (iii '-a). We can 
classify these oscillations into two types, according to the current amplitude. The first 
type is observed in the single plateau (figure 2(i)) and the lower part of the double- 
plateau (figure 2 (iii)) with large amplitude (red). The plateau in this case is located in 
the middle of the NDR region. The second type is observed in the upper part of the 
double-plateau (figures 2(ii) and 2 (iii)) with small amplitude (green). The plateau in 
this case is located at a current approximately three-quarters of the peak current. The 
second type of oscillation contains two Fourier components (figures 3 (ii-a) and (iii-a)). 
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Figure 2. (a) The I-V characteristics for three sizes of the contact regions: (i) 16.950 
nm; (ii) 33.900 nm; (iii) 42.375 nm. The black curve with circles represents the case 
in which the bias is increasing, and the blue curve with the xs represents the case in 
which the bias is decreasing, (b) Time evolution of the current for the value of the 
bias voltage at which current oscillation occurs. The red curves represent the current 
oscillation with large amplitude, and the green curves represent the current oscillation 
with small amplitude. The inset in (ii-b) contains a close-up view of one portion of 
the current. The snapshots of the Wigner distribution at the time points marked with 
purple circles are presented in figures 4 and 5. 



whereas the first type contains just a single component (figures 3 (i-a) and (iii'-a)). 
We find that as the size of the contact regions increases, the frequency of each peak 
decreases. 

As mentioned by Kluksdahl et al. [60] and Zhao et al. [62j, the quantum well formed 
on the emitter side of the effective potential plays an important role in the realization 
of hysteresis and plateau-like structure in the NDR region. The time averaged electron 
densities (black dashed curves) and the effective potentials (black solid curves) calculated 
from the Wigner distribution function in the case of increasing bias are plotted in figures 
3 (i-b)-(iii'-b). For reference, in Appendix A, we present a graph corresponding to figure 
3 (iii'-b) depicting the situation in the case of decreasing bias. A basin-like potential 
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Figure 3. (a) The frequency distribution of the current oscillation that arises in 
the plateau region in the case of increasing bias. The size of contact region and the 
bias voltage are as follows: (i) 16.950 nm and 0.332 V (red curve in figure 2(i-b)); 
(ii) 33.900 nm and 0.280 V (green curve in figure 2(ii-b)); (iii) 42.375 nm and 0.292 
V (green curve in figure 2(iii-b)); (iii') 42.375 nm and 0.302 V (red curve in figure 
2(iii-b)). The insets depict the corresponding transitions. The colored lines represent 
the eigenstates in the emitter basin given in the right graph, while the thick grey 
line represents the continuous energy band, (b) The time-averaged effective potential 
(black solid curve) and time averaged electron density (black dashed curve) for (i)-(iii'). 
The red, green, blue, orange, and purple curves represent the eigenstates in order of 
increasing eigenenergy calculated using the averaged effective potential without the 
heat bath. 



on the emitter side (emitter basin) is observed in the case of increasing bias, while the 
emitter basin does not exist in the case of decreasing bias. In the case of increasing 
bias, when the bias exceeds the peak point of the I-V curve, the kinetic energy of the 
inflowing electron becomes larger than the eigenenergy of the resonant tunneling state. 
As a result, the reflection of the current from the emitter side of the barrier becomes 
large, and the electron distribution function becomes concentrated near the barrier, as 
depicted in figure 3 (i-b) - (iii'-b). When the electron density increases, the effective 
potential decreases. As a result, the emitter basin appears. When the emitter basin 
becomes sufficiently large, there appear resonant tunneling states between the emitter 
basin and the double-barrier well. In such a case, we find current oscillation and a 
plateau of the I-V curve in the NDR region. In the case of decreasing bias, however, 
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because there is no resonant tunneling state between the emitter and double-barrier 
well, the current is much smaller than in the case of the increasing bias. This difference 
causes the hysteresis behavior. 

To elucidate this point more clearly, we solve the Schrodinger equation for the 
emitter basin and the double-barrier well to obtain approximate eigenstates and 
eigenenergies of an electron whose energy is lower than the continuous energy band 
on the emitter side. It should be noted that we employ the time-averaged potential for 
the purpose of the graph, but the effective potential and the corresponding eigenstates 
actually vary in time, because they depend on the electron distribution function, and it 
varies in time. Thus, for example, the identifications of the first (red) and second (green) 
eigenstates and the third (blue) and fourth (orange) eigenstates change in time, often 
becoming degenerate and interchanging. In addition, we ignore the continuous band on 
the collector side and the influence of the heat bath when calculating the eigenstates 
and eigenenergies. For this reason, the calculated eigenenergies are not precise, but we 
find that they are sufficient for determining the cause of the current oscillation, because 
each resonant frequency is rather isolated when estimating the oscillation frequency. 

We find that each oscillation peak in both the large and small oscillation cases can 
be attributed to transitions between eigenstates in the emitter basin, as depicted in 
the insets of figures 3(i-a)-3(iii'-a). As shown in figures 3 (i-b) and 3 (iii'-b), the first 
eigenstate is the tunneling state in the large oscillation case, while the third eigenstate 
in figures 3 (ii-b) and (iii-b) is the tunneling state in the small oscillation case. When we 
compare the profiles of each eigenstate and the electron density distributions, we find 
that the tunneling state and the higher energy eigenstate close to the tunneling state 
are populated in both cases. Thus, we conclude that the current oscillation results from 
transitions between these two states, with the frequency of the oscillation determined by 
the frequency of these transitions. Because the structure of the emitter basin is stable 
with respect to change of the bias voltage, a plateau forms. As the size of the contact 
regions increases, the transition frequencies decrease, because the size of the emitter 
basin increases, while the depth of the basin does not change. This accounts for the 
peak shift seen in figures 3 (i-a) and 3 (iii'-a) and figures 3 (ii-a) and 3 (ni-a). Because 
the emitter basin becomes more stable with respect to change in the bias voltage as the 
contact regions becomes smaller, the plateau becomes large when we decrease the size 
of contact regions. 

Both figures 2 (ii) and (iii) exhibit double plateau-like features in the NDR region, 
but we find current oscillation only in the upper plateau in the case of figure 2 (ii). This 
is because the dissipation in the large oscillation case of figure 2(ii) is not sufficiently 
strong to create significant population of the tunneling states in the case of a small 
basin, for which the resonant frequencies between the tunneling states and adjacent 
states are large. In our previous study [63], we considered the same condition as in 
figure 2 (iii), with a value of 7 half as large (= 12.1 THz). In that case, however, there 
was only upper plateau oscillation. This is because the effective system-bath coupling 
strength, estimated as oc (■j^Uc/{'j'^ + ujI), where Uc is the characteristic frequency of 
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Figure 4. Current oscillation with a large amplitude observed in the lower part of 
the double-plateau structure at the times marked on the red curve in figure 2 (iii-b) 
depicted as snapshots of the Wigner distribution, (a) Current flows into the system 
from the emitter side of the boundary, and then it is scattered (pin — > ~Pim where 
Pin is the momentum of the inflow current) by the emitter side of the barrier almost 
elastically. The scattered electron flows into peak B. (b) Peaks A and B become higher 
and lower in turn, exhibiting motion similar to that of pistons in a two-piston engine. 
(c) The current becomes large whenever peak A becomes large due to the tunneling. 



the system [33], in the previous case is weaker than in the present case. Thus, in that 
case, the ground tunnehng states are not populated from the conduction band through 
dissipation. Note, however, that if the system-bath couphng is too strong, dissipation 
suppresses the current oscillation, as explained in Appendix B. 

One important aspect of the present methodology is that it allows elucidation of the 
dynamical behavior of the system through the time evolution of the Wigner distribution 
function. We have been able to characterize the patterns of the time evolution of the 
Wigner distribution for two types of current oscillations, with large and small oscillation 
amplitude. Here, we describe these two types in detail, using as one reference, the large 
oscillation case in figure 2 (iii-a) and the small oscillation case in figure 2 (ii-a). In 
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figure 4, we display snapshots of tlie Wigner distribution function for tlie case of large 
oscillation at the times marked on the red curve in figure 2 (iii-b). As illustrated in 
figures 3 (i-b) and 3 (iii'-b), the characteristic feature of this type of oscillation is the 
large electron density near the emitter side of the barrier, which is observed as a distinct 
peak separated from the conduction band (at g = 42 nm in figure 3 (iii'-b)). As a result 
of this feature, the effective potential possesses a deep emitter basin next to the barrier. 
The profiles of the eigenstates depicted in figure 3 (iii'-b) indicate that this peak consists 
of the first (red) to fourth (orange) eigenstates. We find that a small peak near the edge 
of the conduction state (at g = 35 nm in figure 3(iii'-b)), which arises from the third 
(blue) and fourth (orange) eigenstates in figure 3 (iii'-b), also plays an important role 
in the current oscillation. In the Wigner distribution plotted in figure 4(a), these two 
peaks are denoted by A and B, respectively. In the situation depicted in figure 4 (a), 
the current flows into the system from the emitter side of the boundary, and then it is 
scattered by the emitter side of the barrier almost elastically, because the kinetic energy 
of the current electron is much higher than that of the tunneling state. The scattered 
current is trapped by the emitter basin and rotates clockwise around the peaks A and B, 
but, due to dissipation, some of it flows into peak B while losing energy. In the eigenstate 
representation given in figure 3(iii'-b), this behavior corresponds to population transfer 
from the fourth (orange) to the third (blue) eigenstate. In figure 4 (b), when the third 
(blue) state decays further to the first (red) and second (green) tunneling states, the 
height of A increases. As a result, peak A becomes higher, while peak B becomes lower. 
In figure 4(c), because peak A is related to the tunneling state, the outflow current 
becomes larger whenever peak A becomes higher. Throughout this process, the peaks 
A and B become higher and lower by turns, in a manner reminiscent of the piston in 
a two-piston engine. As a result of this motion, the current exhibits oscillation. This 
behavior is typical for this large oscillation case. Because there is no current in the case 
of Fig. 4 (a), the oscillation amplitude is large compared to that in the small oscillation 
case. If the dissipation is too strong, however, the piston-like motion is suppressed, and 
there is only steady current, as described in Appendix B. 

In figure 5, we display snapshots of the Wigner distribution for the case of small 
oscillation at the times marked in figure 2 (ii). In figure 5 (a), while current flows into 
the system from the emitter side of the boundary, that part of it with energy larger than 
that of the tunneling state is scattered by the emitter side of the barrier. The remaining 
current, i.e., that whose energy is closer to the energy of the tunneling state, flows to the 
collector side in the form of steady current, through tunneling. In figure 5 (b), it is seen 
how the scattered electron moves back and forth in the emitter basin, while losing energy 
due to dissipation. As a result, the electron flows into peak C, exhibiting tornado-like 
motion. Figure 5 (c) depicts shaking motion of the effective potential that periodically 
accelerates the electron distribution in peak C to the tunneling state. Through this 
effect, current flows to the collector side through the barrier. Due to synchronization 
with this shaking motion, the current is enhanced periodically. This tornado-like motion 
is typical for the small oscillation case. Because there is a large contribution from steady 
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Figure 5. Current oscillation with small amplitude observed in the upper part of the 
double-plateau structure at the times marked on the green curve in figure 2 (ii) depicted 
as snapshots of the Wigner distribution, (a) Current fiows into the system from the 
emitter side of the boundary. Then, a part of the current is scattered (pi„ — >■ —Pm, 
where pin is the momentum of the inflow current) by the emitter side of the barrier. 
The other part of the current flows to the collector side in the form of steady current 
through tunneling, (b) The scattered electron flows in a tornado-like manner to peak 
C in the emitter basin due to dissipation, (c) The shaking motion of the effective 
potential periodically accelerates the component at C to the tunneling state, and the 
current is thus enhanced. 



current, the oscillation amplitude here is smaller than in the large oscillation case. 



5. Conclusions 



In summary, we investigated current oscillations in the plateau structures of the NDR 
region for three sizes of the contact regions with a model that includes damping, 
employing the Caldeira-Leggett Hamiltonian. We found two distinct types of current 
oscillations. The first type is observed in the single plateau and in the lower part of 
the double-plateau structure. It is characterized by a large oscillation amplitude and 
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a single Fourier component. The other type is observed in the upper part of double- 
plateau structure. It is characterized by a small oscillation amplitude and two Fourier 
components. An emitter basin that forms on the emitter side of the effective potential 
plays a key role in creating the current oscillation. Eigenstate analysis indicates the 
that the first type is caused by transitions between the ground tunneling state and 
the adjacent excited state in the emitter basin, while the second type is caused by 
transitions between the intermediate tunneling state and higher states. Because the 
transition frequencies are large in the case of broad emitter basin, there is high frequency 
oscillation in the case of small contact regions for a fixed basin depth. In Wigner space, 
these two types of oscillation are characterized by the two types of motion: two-piston 
engine-like motion and tornado-like motion. Dissipation plays an important role in 
the realization of current oscillation. In order for the ground tunneling state to be 
populated in the case of large oscillation, there must be fairly strong dissipation, whose 
strength is determined by the system-bath coupling and the noise correlation time. If 
the dissipation is too strong, however, the current oscillation vanishes due to damped. 

On the basis of the reduced hierarchy equations of motion (HEOM) approach, our 
investigation was carried out through highly accurate numerical calculations applied 
to a tunneling device system in a non-Markovian environment at finite temperature. 
While previous studies of this kind used the Boltzmann equation [GQl [611 |62], we 
used the Caldeira-Leggett model. Although there is empirical evidence suggesting that 
the Boltzmann equation can be used to treat quantum dissipative dynamics, it is not 
quantum mechanically rigorous, unlike the Caldeira-Leggett model. We thus believe that 
our approach is superior. Although the validity of the Caldeira-Leggett Hamiltonian in 
the description of electron motion is not yet well established, we believe that the present 
results provide insight into the role of quantum mechanical phenomena in the type of 
system studied here. 

The present approach can be used to treat a strong system-bath coupling non- 
perturbatively. In addition, any time- dependent external field can be added while taking 
into account the system-bath quantum coherence through the hierarchy elements. Such 
features are ideal for studying SQUID rings [9l[T0] and quantum ratchet systems [Tl][T2]. 
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Figure Al. (a) The time-averaged effective potential (black solid curve) and time 
averaged electron density (black dashed curve) in the case of decreasing bias. Here, 
the width of the contact region is 42.375 nm, and the bias voltage is 0.302 V. In 
contrast to the cases considered in figures 3 (i-b)-(iii'-b), in the case considered here, 
the emitter basin is very small. The red and green curves represent the first and 
second eigenstates. (b) The steady-state Wigner distribution. The arrow indicates the 
direction of the steady current. 



Appendix A. Effective potential and Wigner distribution in the case of 
decreasing bias 

To understand the origin of hysteresis in the NDR region, we plot the effective potential, 
electron density and energy eigenstates of the emitter basin and double-barrier well in 
the case of decreasing bias in figure Al (a). In this case, the emitter basin is so shallow 
that there is no tunneling state between the emitter basin and the double-barrier well. 
The existence of the peak near the barrier indicates that the second excited state is 
significantly populated. The current arises through the transition from the second 
excited state to the ground quantum state through dissipation. The Wigner distribution 
is plotted in figure Al (b). Due to dissipation, the distribution is in a steady state. The 
peak near the barrier arises because the barrier impedes the flow. The electron density 
then leaks to the collector side without oscillation through tunneling in the form of 
steady current. 



Appendix B. Strong coupling case 

To see the effect of dissipation, we determined the I-V characteristics for the case of 
a stronger system-bath coupling, ( = 120.8 GHz {(^^ = 8.28 ps). The values of the 
other parameters are the same as in figure 2. In contrast to the case considered in figure 
2, in the present case current oscillation is not observed even in the plateau region. 
Also, the upper plateau does not exist. This is because the current oscillation decays 
quickly through the damping. Because the heat bath is strongly coupled to the system, 
the eigenstate picture of the electron system itself, as depicted in figure 3(b), is of 
questionable validity. As a result, the plateau is lost. 
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Figure Bl. The I-V characteristics for the case of strong system-bath coupling, 
C — 120.8 GHz {(^^ — 8.28 ps). The values of the other parameters are the same 
as in figure 2. The black curve with the circles represents the case in which the bias 
is increasing, and the blue curve with the xs represents the case in which the bias 
is decreasing. Comparing this figure with figure 2(a), it is seen that the size of the 
NDR region is smaller and the plateau structure is less pronounced here than in the 
weak-coupling case. Current oscillation is not observed even in the plateau region. 
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